function sort_rate(t_start,t_end,sbin_ms,step_ms,period_ms,neuron_id1,neuron_id2)

%function sort_rate(t_start,t_end,sbin_ms,step_ms,period_ms,neuron_id1,neuron_id2)
% e.g. if you have files up to spikes5000.dat, you have times from
%  0 to 4999, so use something like this: sort_rate(0,4999,100,10,1000,97108,97591);
% t_start time to start the analsis in msec (counting from 0)
% t_end time to end the analysis in msec. (counting from 0)
% neuron_id1 -- starting neuron number to look at (assuming counting from 0)
% neuron_id2 -- end neuron number to plot (assuming counting from 0)
% sbin_ms -- bin size in msec. to calculate spike firing rate
% step_ms -- time step in msec. to shift spike-rate windows (sbin_ms)
% period_ms -- period of behavioral trials in msec.

%t_start=53000;
%t_end=60000;
neuron_id1=neuron_id1+1;
neuron_id2=neuron_id2+1;
 
% read_groups;  % --- load groups.dat, can be commented out for this
% function since celltype information is not used here.


sp=[];
for time=t_start:1000:t_end
    filelabeltime = floor((time+1000)/1000)*1000;
    tmp=load(['spikes' num2str(filelabeltime) '.dat']);
    sp=[sp;tmp];
end;

t_start = t_start + 1;
t_end = t_end + 1;

s = sparse(sp(:,1)+1,sp(:,2)+1,ones(size(sp,1),1),t_end,max(neuron_id2,max(sp(:,2))+1));

%sbin_ms=100;
%step_ms=10;

sec=(t_end-t_start)/1000;

cycles=(t_end-t_start)/period_ms;

for t=t_start:step_ms:t_end-2*sbin_ms

    bt=sum(s([t:t+sbin_ms-1],neuron_id1:neuron_id2),1);
    b((t-t_start+1000)/step_ms+1,:)=full(bt);
end;

t_step=(t_end-2*sbin_ms-t_start+1000)/step_ms-1;


for ii=1:cycles
    bavg(ii,:,:)=b((ii-1)*period_ms/step_ms+1:ii*period_ms/step_ms,:);
end;
bavg2=mean(bavg,1);
bavg2=squeeze(bavg2);
b=b';
bavg2=(bavg2)';
[bmax bmaxind]=max(bavg2');
[btmp bind]=sort(bmaxind);

figure
set(gca,'fontsize',14);
imagesc(bavg2(bind,:));
title(['Neuron # ' num2str(neuron_id1) ' to ' num2str(neuron_id2) ' in one cycle']);
% converting xlabel to msec. scale:
tmp=get(gca,'XTick');
%set(gca,'XTickLabel',tmp*step_ms);
xlabel('Time (msec.)')
colorbar



